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Abstract 

The Wiener-Hopf and Cagniard-de Hoop techniques are employed in order to 
solve a range of transient thermal mixed boundary value problems on the half¬ 
space. The thermal held is determined via a rapidly convergent integral, which 
can be evaluated straightforwardly and quickly on a desktop PC. 
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1 Introduction 

Traditionally, great interest has been shown in determining the disturbances that are 
generated when loads are applied on the surface of a half-space. Lamb |Tj obtained the 
exact solution when an impulsive, concentrated load is applied along a line of the free 
surface of an isotropic linear elastic medium, de Hoop reappraised this problem [2], 
modifying the method originally devised by Cagniard 0. H. leading to the now well- 
known Cagniard-de Hoop (CdH) technique. This method has been used widely since, 
allowing exact solutions to be obtained for a wide range of transient elasticity problems. 
The method can also be useful in order to render solutions into integral forms that are 
rapidly convergent when calculated numerically. 

Transient thermoelastic half-space problems were considered by Danilovskaya [5], Bo- 
ley and Tolins [bj and Achenbach [7] but in these problems the forcing was such that the 
CdH technique was not required. The extension of these problems to inhomogeneous 
media was considered by Baczynski [8] and Parnell {9]. A purely thermal, transient 
problem that employed the CdH method was solved in [10]. The thermoelastic Lamb 
problem was studied by Nayfeh and Nemat-Nasser m who used generalized thermoe¬ 
lasticity in order to retain a finite thermal wave speed, employing the CdH technique to 
determine the solution. 

All of the above problems are of fundamental importance in an array of applications 
where a number of alternative boundary conditions on the surface can arise. What 
appears to be rather lacking in the literature however are studies of transient problems 
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with mixed boundary conditions, where in the context of the thermal problems the 
condition takes the form, e.g. 

dT 

T(0, y, t ) = f(y, t ) for y > 0, -^-(0, 2/> *) = #0/’ *) for 2/ < °> ( L1 ) 

where T(x,y,t ) is the temperature field, f(y,t) and g(y,t ) are two specified functions, 
t is time and with reference to Fig. [Tja: and y are Cartesian coordinates. The half-space 
resides in x > 0 and y runs parallel to the surface, which is defined by x = 0. 

Generally such problems lead to the propagation of a thermal disturbance into the 
half-space. Indeed, such thermal front problems are of importance in a number of appli¬ 
cations including defect sizing m transient thermography ra. solar cell manufacturing 
H and thermal insulation ng. Caflisch and Keller [16], Levine H3 and Satapathy and 
Sahoo [IS] studied front propagation in the thermal context with mixed boundary con¬ 
ditions but in the context of steady problems with applications in quenching. Kozlov et 
ah [19] considered a transient half-space problem with mixed boundary conditions and 
made progress by using cylindrical coordinates due to the special form of the boundary 
condition chosen. 

Mixed boundary conditions are generally difficult to handle even in steady problems 
and analytical or semi-analytical solutions are frequently only possible by the application 
of the Wiener-Hopf method [2Uj. This method exploits the analyticity properties of 
functions in order to yield an explicit or approximate solution in the Fourier transform 
domain. Contour integration then yields the solution in the physical domain. 

Here we shall consider a rather general mixed boundary value problem in the context 
of thermal front propagation and determine solutions using the Wiener-Hopf method and 
Cagniard-de Hoop technique. This problem of mixed boundary conditions of the form 
(II.ID is of particular interest in analyzing the field close to the location of the change in 
boundary condition type, i.e. x — y — 0 in (11.11) . 

We obtain a solution in single integral form by using a deformation of the Laplace 
contour in a similar manner to the Cagniard-de Hoop method. Although it appears that 
we cannot obtain an explicit solution, the solution determined can be evaluated rapidly 
on a desktop PC and therefore it is of great utility due to its general form and its ability 
to circumvent a direct numerical simulation of the problem. Although similar problems, 
involving a discontinuous temperature boundary condition have been considered in the 
building insulation literature, see e.g. Claesson and Hegentoft [21] and Hegentoft and 
Claesson [22] to the authors’ knowledge it does not appear that the solution we provide 
has been written down anywhere in the literature before now. 

In this paper we shall first set out the problem description in section [2] before de¬ 
termining the solution in the transform domain in section [3] In section [4] we describe 
how we deform the Laplace contour onto a steepest descent path, in the manner of the 
Cagniard-de Hoop technique in order to obtain a solution in terms of a single integral 
along the deformed contour path, with an integrand that decays exponentially. In sec¬ 
tion [5] we illustrate the efficacy of the scheme by determining the solution for a number 
of different boundary conditions, with validation provided by finite element solutions. 

2 Problem description 

Assume that the problem under consideration is two-dimensional, being independent of 
z and define the two dimensional half-space domain V = {(x, y) : 0 < x < oo, —oo < 
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y < oo}. We seek solutions to the anisotropic heat equation: 

k fd 2 T „d 2 T\ dT 

V 


pcv 


dx 2 


dy 2 


dt 


( 2 . 1 ) 


where k and k£ are the thermal conductivities (£ > 0) in the x and y directions respec¬ 
tively, cy is the specific heat at constant volume, p is the mass density, t is time and 
T = T(x,y,t ) is the temperature field. We can combine the constants as k = k/(pcv), 
the thermal diffusivity. 

It is convenient to non-dimensionalise the governing equation, using coordinates 
with a “hat” and scale the y coordinate to remove the anisotropy coefficient. Write 
(x, y, T , t) = (x/x*,y/y*, (T - T*)/T*, t/t*), where 


£* = 1 [m], y* = — [m] 


t* = 


[x 


* \ 2 


K 


w, 


( 2 . 2 ) 


and T* is the reference temperature in Kelvin. Upon doing so and “dropping hats” we 
find 

dT 


V T = 


dt 


(2,3) 


where V 2 = d 2 /dx 2 + d 2 /dy 2 . We wish to solve (12.3(1 on the (scaled) domain V with 
boundary dV = dV~ U dT> + as illustrated in Fig. [0 We consider homogeneous initial 
conditions of the form 

T(x>0,y,t = 0) = 0 (2.4) 

and boundary conditions of the form (11.11) but simplify by removing the y dependence, 


i.e. 


dr 

dx 


T = T{x = 0 , 2 / > 0,t > 0) = T 0 / 0 (t), 

dV+ 

dT 

= -z-ix = 0, y < 0, t > 0) = T* 0 g 0 (t), 
&d- dx 


(2.5) 

( 2 . 6 ) 


where T 0 and Tq are real constants and fo{t) and go(t) are piecewise continuous functions 
of time. 

We therefore have a mixed boundary value problem, which in general are not straight¬ 
forward to solve even in the steady context so that the time dependence adds an ad¬ 
ditional element of complexity. Furthermore we allow for the fact that we could have 
a step change at t — 0 on x — 0, leading to a propagating discontinuity front in the 
half-space. 

In order to determine T it is convenient to introduce an alternative problem (for 
convergence issues as will be shown), giving rise to a different temperature distribution 
T, depending on a small parameter e and where T converges to T as e —> 0. This 
problem is described as follows 


dT 

dx 


. m dT 

v ' T = ~m 

T(x > 0, y, t — 0) = 0 
= T(x = 0, y > 0, t > 0) = T 0 f 0 (t)e~ ty 
dT 

= -x -(x = 0,y < 0,t > 0) = T/g 0 (t)e ey . 
dv- dx 


dT>+ 


(2.7) 

( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 


As is easily seen, we recover the solution to our original problem by taking the limit as 
e tends to zero: 


T(x,y,t) = \imT(x,y,t). 

e —>0 


( 2 . 11 ) 
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Figure 1: Domain T> of the problem and its boundaries dT> + and (JD 
on which Dirichlet and Neumann boundary conditions are imposed respectively. 

3 Solution in the transform domains via the Wiener- 
Hopf technique 


Define the Laplace transform in time for any function </>(x, y , t ) by 

COO 

£((j)(x,y,t)) = 4>(x,y,s) = / cj)(x,y,t)e~ st dt 

Jo 


(3,1) 


and hence applying this to the governing scaled equations (I2.7MI2.10H we have 

V 2 f = sT (3.2) 

and boundary conditions become 

T(x = 0, y > 0, s) = f 0 (s)T 0 e~ ey (3.3) 


dT 

— {x = 0, y < 0, s) = g 0 (s)TQe ey . 


(3.4) 


Although s G C the set of complex numbers, for the sake of the analysis to follow we can 
assume it to be real and positive. This allows us to scale the (x, y ) variables to simplify 
the governing equations. The derivation goes through retaining explicit dependence on 
s in the governing equation, but the algebra becomes rather heavy and tedious and 
does not render any greater understanding of the problem; both approaches lead to the 
same result. As such we will rescale x and y in order to eliminate s from the governing 
equation. Thus define 


x 0 = xv s, 

and therefore we obtain 


Vo = yVs 


V 2 f = T 


where Vq = <9 2 /<9xq + d 2 jdy^. The boundary conditions become 

f(x 0 = 0, 2/0 > 0,s) = f 0 (s)T 0 e~ ey 

dT ^ g 0 (s) 

— (x 0 = 0,2/o < 0, s) = — —T n e y . 

OXc) 


(3.5) 

(3.6) 

(3.7) 

(3.8) 
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Next define the Fourier transform in y Q as 


F(T (x 0 ,y 0 ,s)) = Q(x 0 ,a : s) 


f(x 0: y 0 ,s)e iayo dy 0 


( 3 . 9 ) 


and define 0 + and 0 as 

0~(xo ,a,s) = f(x 0l y 0 ,s)e iayo dy 0 and 0 + (x o ,a,s) 

J —CO 


f(x lh y 0 ,s)e“‘ t "‘dy rh 

( 3 . 10 ) 


so that 0 = 0 + 0 + . Applying the Fourier transform to the governing (Laplace 

transformed) equation (13.611 we find that 


0 " = (a 2 + 1)0 


( 3 . 11 ) 


and to the boundary conditions, we find that 


0 + (xo = 0, a, s ) 


<90- 
dx 0 


(x 0 = 0, a, s) 


iT 0 

(a + ie) 


M s ) 

9o( s ) 


(a — ie) y/s 


( 3 . 12 ) 

( 3 . 13 ) 


With reference to Fig. [2J we note that @ + is analytic on fl + = {a e C, Q'(a) > —e}, 
0” is analytic on 12“ = {a G C, Q'(a) < e} and 0 + , 0~ and 0 are analytic on the 
strip S = 12“ D 12 + . The superscript + and — notation thus indicates a.nalyticity in the 
domains 12 + and 12“ respectively. 



Figure 2: Illustrating the domains 12 + , and S on which the functions 0 + , 0 and 0 
are analytic, respectively. 


The solution of (13. lip is 

0(rco, ck, s) = Ai(a, s) exp(—( ck 2 + l) 1 ^ 2 ^) (3-14) 


where the branch of the square root function in the exponent is chosen so that its real 
part satisfies 3ft(a 2 + 1) 1//2 > 0, ensuring that the solution decays as x 0 —> oo. For 
conciseness let us introduce $ + and analytic on IF and 12 _ respectively, as 


$ + (a;, s) = 


<90+ 


dx 


x=0 


and T (a, s) = 


i=0 


( 3 . 15 ) 


Imposing the boundary conditions (13. 12jl - (13. 13p . employing (j3. 15j) and eliminating A 1 
between the two resulting equations, we arrive at 


_ $+ 


K(a)V~ 


K(a)iT 0 ~ _ fFp g 0 (s) 

(a + ie) 0 (a — ie) y/s 


( 3 . 16 ) 


5 





























(3.17) 


The kernel here K(a) = ( a 2 + l) 1 / 2 is easily factorized as K(a) = K~ / K + where 
K + (a) = (ct + z) -1 / 2 and K~(a ) = (a — i ) 1 ^ 2 . 

Multiplying (13.16P by K + we obtain 

-K + (a)$ + = K~(a)V~ + S(a, s) 

where 

K~(a)iTo ~ K + (a)iT^g 0 (s 
S{a,s) = — ——r-/o(s) ■ 


[a + ie) ' ' {a — re) 
= S' - (a, s) + S + (a, s), 




(3.18) 

(3.19) 

(3.20) 


where we have indicated that we wish to determine a sum factorization of the function S. 
One can employ the pole removal method [23] in order to show quite straightforwardly 
that 


3 (a,s) = iT 0 f 0 (s)L 1 + 


v9o( s ) r + 

0 x/i 2 ’ 

(3.21) 

,f9o(s) 

0 V~s 2 ’ 

(3.22) 


where 


L~ x = 


Lt = 


a 


i ) 1 ' 2 


C- 


a + ie 

(ck + i ) -1 / 2 — 
a — ie 


c+ 


r+ _ c - 

1 a + ie” 

T -- C + 

2 a — ie 


and 


c_ = 


— [— i(e + 1)] 


1/2 


c+ = [i(e + 1)] 1/2 . 


Referring to (I3.16j) we can therefore define a function ^(a) such that 

E(a) = 


-(K+&+ + S+) 
-(K+$ + + S + ) = + S" 

(k~v- + s~) 


on 0 + \iS 
on 5 
on Q~\S, 


(3.23) 

(3.24) 

(3.25) 

(3.26) 


and therefore is analytic on the whole a-complex plane. It can be shown that as |a| —y oo, 
E(a ) = 0{a^ 1 ^ 2 ) when a E and E(a) = 0(a 1//2 ) when a E £l~. As such E(a) = o(a) 
as |a | —y oo. This together with the analyticity of E(a) and the extended Liouville 
theorem (see for example [20]), implies that E(a) is constant. However, we also know 
that E(a) —> 0 as |a| —> oo and a E Q + . We can then conclude that E(a) = 0 
everywhere. Given this, we therefore have from (13.261) 


$+ = 


S+ 

W+ 


and T = 


S~ 


K~ 


(3.27) 


From the original expressions for A 1 determined from the boundary conditions we can 
show that 


Ai(a,s) = 


iTn 


fo(s) + 


(a + ie) 

iT 0 c-fo(s) 

(a + ie)(a — i) 1 / 2 ' (a — ie)(a — i) l / 2 ^/s 


+ 


iT{,c + ~g 0 (s) 


(3.28) 

(3.29) 
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where (13.17j) . 03.23p . (I3.24p and (I3.22p have all been used. Referring to (I3.14p . the 
solution in transform space is therefore 


0(x o ,a, s) = 


iT 0 ff 0 {s ) 


+ 


iToggo(s) 


(a + ie)(a — i) 1 / 2 (a — ie)(a — i) 1 / 2 ^/!. 


ex' 


p(-(a 2 + l) 1 / 2 x 0 ). 


(3.30) 


Formally inverting the Fourier transform and using (13.5p gives the Laplace trans¬ 
formed solution as 


\ iT 0 fo{s) T t \ i OZffo(s) r , \ 

T{X, y, a) = —-- h(x, y, s) + — —^/ 2 (x, y, s), 


2tt 


where 


h = c_ 


roo g—[(a 2 +l ) 1 / 2 x+iay\-</s 

Loo ( a + *e)(a — i) 1 / 2 


da and J 2 = c + 


27Ta/s 

/ oo g—[(« 2 +l) 1 / 2 x+iQ!j/]-v/5 

-oo (at ~ ie)(a - i) 1 / 2 


(3.31) 


da (3.32) 


4 Semi-analytical inversion via a Cagniard-de Hoop 
approach 

Motivated by the Cagniard-de Hoop technique, let us introduce polar coordinates r and 9 
related to x and y in the usual manner, i.e. x = r cos 8,y = r sin 9 where 6 G [— 7t/2, 7t/2] 
and introduce the parameter /3 via the expression 

/3r = (a 2 + l) 1,/2 x + iay (4.1) 


so that 


P = (a 2 + 1) 1//2 cos 6* + ia sinR (4.2) 

Inverting for a we therefore determine the two paths B + and B_ in the right and left 
halves of the a-plane 

a± = —ip sin 9 ± \JP 2 — 1 cos 9. (4.3) 

In the a-plane, with P e [1, oo) these paths start at a = —i sin 9 and move off to 
infinity either in the upper half-plane (9 e (—tt/ 2, 0), see Fig. [3]) or the lower half-plane 
(9 e (0,7 t/ 2), see Fig. H]). 

Since B ± are steepest descent paths for the integrals, the idea is to deform the 
integrals (I3.32p from the real line onto these to aid convergence. In classical Cagniard 
de Hoop problems this frequently permits one to write the a integral in the form of a 
Laplace transform of a function that is independent of s and thus we can determine 
the inverse transform immediately, thus rendering explicit solutions. Here we are not 
so fortunate, the function will not be independent of s but nevertheless we are able to 
make significant progress due to the fact that the inverse Laplace transform integral can 
be determined analytically in many important cases as we shall see shortly. This leaves 
the solution in a single integral form that is rapidly convergent. 

At this point note that 

lim c_ = (e"^/ 2 ) 1 / 2 = e _i7r/4 = c and lim c + = ( e i7r/2 )- 1/2 = e" i7r/4 = c (4.4) 

e—s-0 e—>0 v 7 v 

and let us consider the case of negative and positive 9 separately. 
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Figure 3: Illustrating the location of the paths B* 2 when 9 e [—7r/2, 0). The poles at 
a = ±ie and branch point at a — i are also depicted where the branch cut passes to 
infinity vertically along the imaginary axis. 


4.1 The case of 9 £ [— 7r/2,0) 

4.1.1 Evaluation of I\ 

Referring to 03.321) and Fig. [31 we see that the integrand of I\ has a pole at a = —ie and 
a branch point at a — i. Deforming the contour from the real axis to the B ± contour in 
the a-plane we do not cross any singularities and hence we find that 


= c_ 


[(o 2 +l) 1 / 2 x+iay] x /5 


da — c_ 


' b- 


g-[(a 2 +l) 1 / 2 a:+ia 3 /] x /i 

(a + ie)(a — i) 1 / 2 


da 


l B + (a + ie)(a - i ) 1 / 2 

= It - J i 

On B±, we have a = a± = —i/3 sin 6 ± y//3 2 — 1 cos 6 and so we have 

da = (—i sin 0 ± (3 cos 9(/3 2 — l) -1 / 2 )d/h 

Noting that f3 e [1, oo) is a parametrization of the paths B ± we can rewrite if as follows, 
taking the limit as e —)■ 0, 

3 H sin 0±/3 cos^CS 2 -!)- 1 / 2 ) 


(4.5) 

(4.6) 


lim if = c 


e—>0 


a±(a± — i ) 1 / 2 


~^df3 


Therefore from (14. 5 p we determine the form 


lim It = c 

e —>0 


T(p,e)e~ flr ^ d/3, 


(4.7) 


(4.8) 


where 


J-((3, 9) = — ism9 


+ 


P 


a + (a+— i) 1 / 2 a_(a_ —i) 1 / 2 

1 1 


cos 9 


+ 


a + (a_|z) 1 / 2 a_(a_—z) 1 / 2 


(/3 2 - 1)V2 

As an aside, we note by deforming into the lower half-plane that 

0 1 


J-((3, 9) d/3 = lim 


e->°(a + ie)(a — i) 1 / 2 
2in 


da 


(4.9) 


(4.10) 


c 


(4.11) 
















4.1.2 Evaluation of /? 


Referring to (13.321) and Fig. [3j we see that the integrand of R has a simple pole at 
a = ie and a branch point at a — i. Deforming the contour from the real axis to the 
B ± contour in the a-plane we cross the simple pole and pick up its residue. Accounting 
for this contribution we find that 


h = c+ 


3 — [(ct 2 + l) 1 / 2 x+ia:j/]- v /s 


B + (a — ie)(a — i) 1 / 2 
= If — / 2 + R 2 


da — c + 


r g-[(« 2 +l) 1 ' /2 a;+iay]yi 

l B - (a — ie)(a — i) 1 / 2 


da T R 2 


where 


R 2 = 2n ic 4 


3 —(((*e) 2 +1) 1 / 2 cos 9-\-i(ie) smd)ry r s 


lie — l 


0 1/2 


and in particular, 


lim Ro = 2me~ x ^ s . 


e—>0 


(4.12) 

(4.13) 

(4.14) 


Noting that P 6 [1, 00 ) is a parametrization of the paths B ± we can rewrite if as follows, 
taking the limit as e —)■ 0, 

lim It = c (4 . 15 ) 

^0 J 1 a ±( a ± — i ) 1 / 2 

Therefore from (I4.12|) and (14.14j) we determine the form 

roo 

lim I 2 = c R(p, 9)e~^ s dp + 2me~ xV ~ s . (4.16) 


e-s>0 


4.1.3 An expression for T(r, 0, t) 

We show in Appendix [A] that 

cF{p,0) = ^Q{M), 

where Q is a real-valued function. As such, using this together with (13. 31 j) . 
and (14. 17j) we hnd the following expression for T = lim e ^. 0 T 

T(r,9,s) = lim R (x, y, s) + - lim I 2 (x, y , s) 


2n mo 


2ny/s e->o 


V2n 

over 
T(r,9,t) = 


Tofo(s) + R 


Ms) 


g(P, 9)e~^ s dp - T ^M) e -xV~s 


Finally we recover T by taking the inverse Laplace Transform, 
To 


M 


7T J 1 


G(P,9)R(r,P,t) dP + 


T' 
2 0 


(4.17) 

gT6} 


V2i\ J 1 
Tp7^(r, cos 6>, t). 


g(P,9)T 2 (r,P,t) dP 


(4.18) 


where 




/» < 7+200 

<r—ZOO 
/* < 7+200 


fo(s)e-P r M st ds, 

9o{s) c - 0 r ^ c s t 

M 


(4.19) 

(4.20) 

(4.21) 


and where as usual a £ M is chosen here such that all singularities of the integrands are 
to the left of the line s = er. 
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4.2 The case of 9 G (0,7r/2] 

This case follows entirely analogously to the negative 8 scenario, the difference here 
being that the paths B± reside in the lower half of the complex ct-plane and as such the 
pole that leads to a contribution to the integral is at a = —ie. We find that 

rji pOO rpf poo 

T(r,6,t) = J 5(/S,0)7;(r,/S,i)d/3+-J-y 5(/3,0)T 2 (r,/3,()d/3 

+ T 0 T\_(r, cos 8, t). (4.22) 



Figure 4: Illustrating the location of the paths B ± when 8 £ (0,7r/2]. The poles at 
a = ±ie and branch point at a — i are also depicted where the branch cut passes to 
infinity vertically along the imaginary axis. 

As an aside, for 8 £ (0, 7t/ 2], analogously to the derivation of 

TO 3,8) d^ = 0 

Therefore using, (14. lip . (14.171) and (14.231) we have 

-F p g(p. e) = 1 - H (9) 

where H{8) is the Heaviside step function. 

4.3 A summary of the solution 

Combining the results from sections 14.11 and 14.21 we can write down the solution for all 
values of 8 , 

rji pOO rpf pOO 

T(r,e,t) = - 2 = / g((i,e)T l (r,Mdp + -^= j 

+ H($)T 0 Ti (r, cos 0, t) - (1 - H(0))r„T 2 (r, cosff, (). (4.25) 

Both integrals in (I4.25p and the additional term have a discontinuity at 8 = 0 and as 
such this form is not particularly “clean”. We are able to improve upon this form, using 



(14.lip we have 

(4.23) 

(4.24) 
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(I4.24p to generate 


rji POO 

T(r, 0, t) = -2= j S(ft 0){%(r, ft t) - ft(r, cos 0, t)}d/} 

rj~\f POO 

H- 7 =/ G(P,9){T 2 {r,P,t) - T 2 (r, cos6,t)}d/3 + T 0 Ti(r, cos9,t). (4.26) 

7Tv2 J 1 

Each term in this expression is now continuous across 9 = 0. We shall discuss this aspect 
further in the context of specific examples in the next section. 


5 Some specific boundary conditions 

5.1 Perfect insulator on y < 0 

Let us now assume that Tg = 0 so that we have a perfect insulator on y < 0. We shall 
consider a variety of temperature profiles for y > 0. The solution is therefore obtained 
by setting Tg = 0 in (j4.25j) or equivalently fl4.26[) . As such only 71 enters the analysis. 


5.1.1 Step temperature change 

Take the simplest form, f 0 (t) = 1 so that /o(s) = 1/s and we need to determine 71 
defined in (I4.20|) . It transpires that it is convenient to differentiate the expression for 7i 
with respect to t, which enables the inverse Laplace integral to be evaluated analytically 
in this case 


dTi 


dt 




1 

27 n 


PCT+ZOO 


e -prV~s e st 


ds 


r @ r -r 2 p 2 /m 
2^/2' 


(5.1) 


We then integrate (definitely) with respect to t with a lower limit of t = 0, finding that 


Ti(r,P,t) 



Appealing to (I4.25|) we have as our solution 


(5.2) 


T(r, 0, t) 


Tp 

VZlT 



G(P,6) Erfc 



d/3 + 7f(6 l )TgErfc 


r cos 9\ 
2ftf ) 


(5.3) 


noting that G(/3,9) is the real function defined in (IA.12D . Each term of (15.31) is easily 
computed numerically, and, in Fig. [5] we plot the resulting temperature profile on the 
horizontal axis against y running vertically, at time t = 0.02 for two values of x, x = 0.05 
and x = 0.2. The circles are results taken from a finite element solution of the same 
problem in COMSOL and provide validation of the present semi-analytical scheme. 

Note that the temperature profile is continuous across the x-axis, i.e. 9 = 0. Rather 
interestingly, both terms on the right hand side of (15.31) are discontinuous across 9 = 0, 
as is shown in Fig.[6]but the two discontinuities compensate exactly to yield a continuous 
temperature profile. If instead of the form (14.251) . we use (14.261) . the solution is written 
as 


T(r, 9, t ) 


To 

iry/2 


/ OO 

G(P,9) 



— Erfc 


f r cos (9 )^ 

\~WTJ 


d/3 


+ Tq Erfc 


r cos (9 )\ 

2 Vi ) ' 


(5.4) 
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T(x,y) 


Figure 5: Thermal field T(x,y,t) at x = 0.05 (red) and x = 0.2 (blue) when t = 0.02 
for the case of a perfect insulator on y < 0 and a step change temperature increase at 
t — 0 on y > 0. Circles are predictions from solutions to the same problem using finite 
elements methods in COMSOL, which provides validation of the present semi-analytical 
method. 



Figure 6: Plot of the first (left) and second (right) terms on the right hand side of (j5.3jl 
for r = 0.05, t = 0.02 and 6 e [—7 t/2, vt/2]. 

This expression is also straightforward to evaluate numerically and (obviously) gives the 
same results as those presented in Fig. [5j but this time, as shown in Fig. [7J both terms 
on the right hand side of (15.4(1 are continuous across 6 = 0. 

Finally, in Fig. [SJ we plot the two dimensional temperature contour profile on the 
(x, y) plane at t = 0.02 illustrating how the distribution spreads out from the upper 
half-plane. 

5.2 Imperfect insulator on y < o 

Let us now consider the case when Tg ^ 0 and let us take /(f) = g(t) = 1 so that 
/o(s) = go(s) = 1/s and now both 7i and 72 play a role. Once again it is convenient to 
differentiate with respect to f in order to evaluate the inverse Laplace transforms and 
subsequently integrating these expressions definitely with respect to f with a lower limit 
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Figure 7: Plot of the first (left) and second (right) terms on the right hand side of (15.4p 
for r = 0.05, t = 0.02 and 0 G [—7r/2, tt/2]. 

0.4 


0.2 


0.0 


- 0.2 


- 0.4 

0.0 0.1 0.2 0.3 0.4 

Figure 8: Contour plot of the temperature profile when t = 0.02, Tq = 1,Tq = 0, i.e. a 
perfect insulator on y < 0 and step change in temperature on y > 0. We see how the 
thermal field propagates into y < 0. 

of t = 0 yields 



T\ (■ r , p,t)= Erfc (j^=) , ( 5 - 5 ) 

T 2 (r, p, t) = 2^-e-^ 2 /^ - r ?Erfc . (5.6) 

Either of the expressions (14.251) or (I4.26j) then recover the temperature prohle. Both 
formulations are easily computed numerically and give rise to a continuous temperature 
prohle, as seen in Fig.[2]where the prohle is plotted at t — 0.02 for x = 0.05 and x = 0.2. 
The contour plot of the thermal held at t — 0.02 is given in Fig. (TUJ 
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Figure 9: Thermal field T(x, y, t) at x — 0.05 (red) and x = 0.2 (blue) when t = 0.02 for 
the case of an imperfect insulator on y < 0 and a step change temperature increase at 
t — 0 on y > 0. Circles are predictions from solutions to the same problem using finite 
elements methods in COMSOL, which provides validation of the present semi-analytical 
method. 



Figure 10: Contour plot of the temperature profile when t = 0.02, T 0 = 1, TJ) = 1, i.e. an 
imperfect insulator on y < 0 and step change in temperature on y > 0. We see how the 
thermal field propagates into y < 0. 

5.3 Continuous ramp up and down superposed on a step change 

in y > 0 

Thus far we have considered only cases where / 0 and g 0 are constant, accommodating 
for the step change at t — 0 of course. Let us now consider the case when these functions 
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can be unsteady and in particular when f 0 (t) is a step and superposed general ramp up 
and down profile given by 

f 0 (t) = H{t) + {t- a)H(t - a) + 2(6 - t)H(t - 6) + (t - (26 - a))H(t - (26 - a)) 

(5.7) 

and illustrated in Fig. fill 



Figure 11: Plot of /o(t) as defined by (15.71) . A ramp up and down is superposed on a 
step change unit temperature profile. 


The Laplace transform of / 0 (i) is 

r , , 1 , e~ 

fo(s) — - H- 


p—bs g—(2 b—a)s 

— 2—— -h - 77- 


= /o (1) (s) + fo 2 \ s ) + /o 3) (s) + /o 4) (s) 

and the resulting expression for 7i is 


(5.8) 


r*CT+ZOO 


ri (r .ft*) = 5- 

— -tP)/ 


(/o (i) (s) + /o (2) ( s ) + fo 3 \ s ) + /d 4) ( s ))e 


m. 


F(3) 


'( 4 )/ 


= 7?>, p, t ) + 7? 2) (r, A t) + rr(r, t) + mr, /3, t ). 


<3), 


'(4)/ 


(5.9) 


The case of T( (1) has already been dealt with in Section 15.1.11 and is thus given by (15.21) . 
Moreover, T} 2 \r, /3, t), T^ 3 \r, /3,t) and 7^ <4 ^(r, /3, t) have a very similar structure, so 

we only need to consider only the case of 7p (r, /3,t) in detail. In fact, we have 


7? 2) M,t) = I f^\s)e~^e st ds 


f*(T+ZOO 


( 2 ), 


2in 

1 

2in 


rcr+io o g —as 


e~ 0r ^e st ds 


(5.10) 


Differentiating T L (2) twice with respect to time, we obtain 


d 2 Ti 2) 

<9t 2 


(r, /3,6) = H(t — a) 


r/3 




2\fn t(Z — a) 3 / 2 


g 4(i — a) 


(5.11) 
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T(x,y) 


which, following the same reasoning as in Section 15.1.11 implies that 


or. 


( 2 ) 


dt 


■(r, j3, t ) = H(t — a) Erfc 


r/3 


2\Jt — a 

and by dehnite integration with respect to time, we obtain 

r/3 


(5.12) 


T, 2 \r ,/?, t) = H(t — a) <( ( t — a + 


r 2 (3 2 


Erfc 


2 y/t — a 


r(3\/t — a r 2 g 2 

—e 44 


7T 


(5.13) 


The functions r, 3 \r, /3 : t) and r, A \r, /3, t) are obtained in the same manner. Using 
these expressions in the general formulation (14.251) or (I4.26p enables the solution to be 
computed rather rapidly. In Fig. [12] we plot the resulting thermal field at the locations 
(x, y) = (0.1, 0.03) and (0.2, 0.03) in the cases when T 0 = 1 with Tq = 0 (left) and T 0 = 1 
with Tg = 1 (right). We also plot the associated solutions determined previously where 
no ramp up and down is present in the boundary condition. 




Figure 12: Illustrating the effect of the ramp up and down on the evolution of the 
temperature profile at two different locations (x,y) = (0.1,0.03) (red) and (x, y) = 
(0.2,0.03) (blue) for T 0 = 1 and Tq = 0 (left) or T 0 = 1 and Tg = 1 (right). The dashed 
lines correspond to the profile without the unsteady ramp, and the plain lines correspond 
to the profile with the unsteady ramp. 


6 Concluding remarks 

Employing the Wiener-Hopf technique and a Cagniard-de Hoop-type integral, a rapidly 
convergent integral expression has been determined for a class of transient thermal mixed 
boundary value problems. The integral is easily computable on a standard desktop PC 
for a wide range of transient boundary forcings of interest. Here we illustrated the 
computation for a number of cases, including step-changes in temperature and ramp 
up and down boundary profiles. Such quasi-analytical expressions are of great utility 
in order to speed up computations and enable asymptotic analysis close to locations 
of interest. Future work could include extensions to full elastodynamics and coupled 
thermoelasticity. In these cases matrix Wiener-Hopf problems will result in general. 
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A Simplification of the function F(/3,6) 

From (14.9p we have 

1 da+ 1 dcm 




a + (a + — i) 1 / 2 d[3 «_(«_ — i) 1 / 2 dj3 


and 


da± . . (3 

—— = — l Sill 9 ± — , COS 9. 

dp 


Therefore 


F{P,0) = 


a + (a + — i) 1 / 2 




-i sin 9 H- , ^ cos 9 


(A.l) 


(A.2) 


+ 


1 


i sin 9 4- -=t =cos9 I (A.3) 


«-(«- - i) 1 / 2 \ ^//3 2 - 1 


and simplifying further we obtain 


F(P,9) = 


ct_a_|_(a + — i) l R{a_ — i) 1 / 2 


(a + - i) l/2 ((3 


i cos 9 sin 9 


(a_ - i )^ 2 I /3 + 


i cos 9 sin 9 


Next with reference to Fig. [IS] define R and 0± such that 

a + -i = Re^+, a.-i = Rem¬ 

and since 0_ = —tt — 0 + we can write 

(a+ - i) 1/2 = y/Re^+Z 2 , (a_ - i) 1/2 = VRe~^~ /2 = -iVRe~^ 2 . 

Therefore it is possible to show that 

1 


(A.4) 


(A.5) 


(A.6) 


n^) = 


VR{ 1 + i) f3 (cos (0+/2) + sin (0 + /2)) - 


a_a+(a + — 0 1 / 2 (q!_ — i) 1 / 2 

cos0sin6 l ( / 3 2 — 1) _1 Z 2 (cos (0+/2) — sin (0 + /2)) 


Finally, noting that 


and 


a + a- = —(/3 — cos 9) 


(A.7) 


(A.8) 


(a+ — i) 1 Z 2 (a_ — i) 1 / 2 = —i\{3 + sin0| 


(A.9) 
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so that 


a_a + (a + — z) 1//2 (a_ — i ) 1 ^ 2 = i(/3 2 — cos 2 6)\/3 + sin0| 
we have from (IA. 7[) 

jf(/3,») = — g(p,e) 

1C 


(A-10) 


(A.ll) 


where 

SGM) = 


(3 (cos ("0+/2) + sin (ij> + /2)) - 


(/3 2 — cos 2 0)|/3 + sin#! 1 / 2 L 

cos0sin0(/3 2 — l) -1 / 2 (cos (^ + /2) — sin (^ + /2)) 


(A-12) 


is a real valued function. 




Figure 13: Diagrammatic description of the angles 


B Variational formulation used for finite element so¬ 
lution 

In order to compute the transient solution T(x,y,t) of the problem (j2.4p - fl2.6p in the 
semi-infinite domain T> using finite element method, we define a rectangular domain 
Vq e M 2 : V o = [0,a] x [—6/2, 6/2] with boundary dD 0 = dD q" U ODq U dD g° where 
8Dq = {x = 0,—6/2 < y < 0} and 3Dq = {x = 0,0 < y < 6/2}. The parameters 
a, 6 which define the size of T> 0 are chosen sufficiently large such that the temperature 
field on dD “ is not influenced by the perturbation due to the discontinuous boundary 
condition on 8Dq and ODq in the time interval under consideration. As a consequence, 
the symmetric boundary condition may be imposed on dD^. The problem (12.4l) - (i2.6H 
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may be rewritten as: 


f- v ' 2r =°- 

in dD 0 

(B.l) 

T = T 0 f 0 (t), 

on dVf 

(B.2) 

VT • n = -T* 0 g 0 (t), 

on dDf 

(B.3) 

VT ■ n — 0, 

on dD “ 

(B.4) 

e outward pointing normal vector to 8Dq. 



Let 5T be the trial function of T, the weak formulation of this problem reads: For 
any test function 5T G V(V 0 ), find T G V(V 0 ) such that 



5T ^dv + 
at 


V(<5T) • VTdv 


i v 0 


'dV~ 


(ST)T^g 0 ds, 


(B.5) 


and T = T 0 f 0 on dVf, (B.6) 

where V(V 0 ) = { f(x,y ) G H\V 0 )}, V(V 0 ) = { f(x,y ) G H\V 0 )- f{x,y) = 0 ,(x,y) G 
ST>+}. 

This weak formulation has been implemented in the hnite element software COMSOL 
Multiphysics [23] • For simulations in the time interval t = [0, 0.02] as presented in Section 
5, a rectangular domain ( x,y ) G [0,1] x [—1,1] has been used. The element size (h e ) 
and the time step (At) needed for the discretization are respectively h e = 0.01 and 
At = 10~ 4 . 
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